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ABSTRACT 

We combine a semi-analytic model of galaxy formation with simple analytic 
recipes describing the absorption and re-emission of starlight by dust in the inter- 
stellar medium of galaxies. We use the resulting models to predict galaxy counts and 
luminosity functions from the far-ultraviolet to the sub-mm, from redshift five to the 
present, and compare with an extensive compilation of observations. We find that 
in order to reproduce the rest-UV and optical luminosity functions at high redshift, 
we must assume an evolving normalization in the dust-to-metal ratio, implying that 
galaxies of a given bolometric luminosity (or metal column density) must be less extin- 
guished than their local counterparts. In our best-fit model, we find remarkably good 
agreement with observations from rest ^ 1500 A to ~ 250 /im. At longer wavelengths, 
most dramatically in the sub-mm, our models underpredict the number of bright galax- 
ies by a large factor. The models reproduce the observed total IR luminosity function 
fairly well. We show the results of varying several ingredients of the models, including 
various aspects of the dust attenuation recipe, the dust emission templates, and the 
cosmology. We use our models to predict the integrated Extragalactic Background 
Light (EBL), and compare with an observationally-motivated EBL model and with 
other available observational constraints. 
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1 INTRODUCTION 

New observational facilities have greatly extended the range 
of the electromagnetic spectrum over which emission from 
galaxies can be measured, while simultaneously expanding 
the range of cosmic history that can be probed. For ex- 
ample, in recent years, survey observations over significant 
fractions of the sky have been carried out in the Far and 
Near Ultra-violet by GALEX (Galaxy Evolution Explorer), 
in the optical by the Sloan Digital Sky Survey (SDSS), in the 
Near- infrared (IR) by 2MASS (Two Micron All Sky Survey), 
in the Near-IR and mid-IR by the Spitzer Space Telescope, 
and most recently, in the far-IR by the Herschel Telescope. 
In addition, multiple ~ 0.5 — 1 sq. degree sized fields have 
now been deeply imaged in the X-ray with Chandra, in the 

* E-mail: somcrville@stsci.cdu 



optical with the Advanced Camera for Surveys (ACS) on 
the Hubble Space Telescope (HST) as well as ground-based 
facilities, and in the near-IR with UKIRT and Spitzer, al- 
lowing large samples of high redshift galaxies to be identi- 
fied and studied. The recently installed Wide Field Camera 
3 (WFC3) on HST is in the process of observing many of 
these fields at high resolution in the Near IR. Crucial to the 
extraction of physical quantities and scientific insight from 
these deep surveys has also been the availability of accu- 
rate redshift information for large numbers of galaxies from 
multi-wavelength medium-band surveys (COMBO-17, COS- 
MOS, MUSYC, NEWFIRM) and multi-object spectroscopy 
(DEEP, VIMOS). 

However, perhaps some of the most surprising and 
poorly understood observational results of the past two 
decades have come from long-wavelength observations in 
the mid- to far-IR. The IRAS satellite revealed that 
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^ 30% of the bolometric luminosity of nearby galaxies, 
mostly normal spirals, is r eprocessed by dust in the IR 
l|Soifer fc Neugebauerl Il99ll '). and discovered a population 
of heavily obscured luminous a nd ultra-luminous infrare d 
galaxies (LIRGs and ULIRGS: [Sanders fc Mirabell ll996D . 
IRAS already provided hints of the very strong evolution 
of this IR-bright population, the number density of which 
seems to have been enormously larger in th e past. This 
was c o nfirm ed and quantified first by IS O (e.g. Elbaz et alj 
19991 . 2002"), and then by SCUBA Isinail et al.1 ll997F 
Hughes ct al. 199 IChapman et al]|2005l) an d Spitzer (e.g. 
Le Floc'h eral]|2005l : iBabbedge etalll2006l ). The physical 



interpretation of these high redshift LIRGS and ULIRGS, 
however, has until recently been hampered by the fact that 
in many cases observations existed only in the mid-IR or 
sub-mm, relatively far from the ~ 100 /im peak of the dust 
emission (e.g. 15 fim in the case of ISO, 24 ^m for Spitzer, 
and 450 and 850 nm in the case of SCUBA). This situation 
should improve greatly in the next few years, as observations 
bracketing the 100 fim peak are taken by the PACS (57 to 
210 /im) and SPIRE (250, 350, 500 /im) instruments on the 
recently launched Herschel telescope. 

This rainbow of observations presents a challenge to 
theoretical models of galaxy formation. To date, most stud- 
ies have focussed on making predictions for rest-optical or 
intrinsic (e.g. stellar mass, star formation rate) properties 
of galaxies, mainly because of the difficulty of modeling the 
absorption and emission of light by dust in the interstel- 
lar medium of galaxies. However, observational estimates of 
intrinsic physical properties from multi-wa velength photom - 
etry suffer from poorly constrained biases l|Lee et al.ll2009l ). 
and quantities such as the star-formation rate (SFR) can dif- 
fer greatly depending on which observational tracer is used 
to estimate them. Therefore, in order to interpret the zoo of 
galaxies selected at different wavelengths (e.g. LBGs, EROs, 
DRGs, DOGs, BzKs, BM/BX, SMGs, etc.), it is important 
to develop models that can accurately predict observable 
quantities over the full range of wavelengths probed by mod- 
ern panchromatic surveys. 

Important theoretical advances have been made in the 
past few years, as well, with the development of radiative 
transfer (RT) codes that, coupled with a model of the dis- 
tribution of stars, gas, and dust in a galaxy, can produce 
detailed panchromatic prediction s of the galaxy's appear- 
ance and photometric properti es l|Silva et al.lll998l : IJonssonI 
I2OO6I : Ijonsson et aklfiood . I2OI0I ) . One approach is to use ide- 
alized galaxy models (e.g. sphe roid plus disc), cou pled with 
a radiative transfer code, as in ISilva et al.1 (| 19981 ). Another 
is to use hydrodynamic simulations to provide more detailed 
spatial and morphological information, as in lJonsson et al.l 
l|2006l ). 

While powerful, these tools are still computationally 
quite expensive. Producing panchromatic predictions for 
statistical samples of galaxies in a cosmological context re- 
mains beyond reach, without some shortcuts. Semi-analytic 
models (SAMs) of galaxy formation, which apply sim- 
ple but physically motivated recipes for the physical pro- 
cesses that shape galaxy formation, within the framework 
of structure formation predicted by ACDM (A Cold Dark 
Matter), can provide predictions of bulk galaxy proper- 
ties (such as star formation and chemical enrichment his- 
tory, radial size, total stellar mass or luminosity, ratio of 



spheroid to disc, etc.) for very large numbers of galaxies. 
SAMs have been shown to reproduce many observed prop- 
erties of galaxies (e.g. Kauffm ann e t all 19931: Cole et akl 
1994 ISomerville fc Primackl [l999; K auffmann et al.l Il998l : 



Cole et al. I I2OOOI : ISomerville et al.ll200ll ). and to agree rea- 



sonably well with the results of numerical hydrodynamic 
simulations in their predictions for basic quantities such 
as t he rate of accretion of cold gas and gal axy merg- 
ers (jYoshida et al.1 I2OO2I : ICattaneo et al.1 |2007| ). In par- 
ticular, recent models that include "radio mode" feed- 
back from Active Galactic Nuclei (AGN) reproduce quite 
well the global properties of massi ve galaxies over a 
broad range of c o smic history (e.g. ICroton et al.l 20061 : 



Bower et al.l 



Monaco et al 



20061: iMenci et al.1 I2OO6I: iKang et al.l I2OO6I : 
' 2OO7I : ISomerville et al.ll2008ah . although re- 



producing the_^ropertiesofJow-m^^ a 
challenge l|Fontanot et al.ll2009al : TGuo et al.ll2010l ). 

Using a set of recipes for gas accretion and cooling, 
merging, star formation, stellar feedback, chemical enrich- 
ment, and optionally black hole growth and AGN feedback, 
a SAM outputs a distribution of ages and metallicities for 
all the stars within the spheroid and disc components of 
a galaxy (more details are given in Section [5}. This infor- 
mation is co nvolved with "simple stel lar population" (SSP) 
models (e.g. iBruzual fc Chariot! |2003| ). which specify, for a 
given stellar Initial Mass Function (IMF), the luminosity 
as a function of wavelength for a stellar population of a 
given age and metallicity, in order to predict the unattenu- 
ated Spectral Energy Distribution (SED) of starlight in the 
galaxy. The predictions of stellar population models from 
different groups have largely converged in their predictions, 
particularly in the UV and optical, making this component 
of the modelling relatively robuslQ. 

For the more difficult step of predicting how this 
starlight is absorbed and re-radiated by dust, one possible 
approach is to couple the prediction s of a SAM d i rectly 
with a radiative transfer code (e.g. I Cranato et al] |200G| : 
iBaugh et a"Ll l2005: Font anot et al.ll2007l ). Because SAMs are 
not able to track the detailed internal structure or morphol- 
ogy of galaxies, this requires the assumption of an idealized 
geometry such as a spheroid plus disc (where the sizes and 
masses of the components are specified by the SAM). How- 
ever, even this approach is prohibitively expensive for large 
numbers of galaxies, and also has the disadvantage that the 
simplified geometries may not be representative of the di- 
versity of galaxy types, particularly for LIRG or ULIRG-like 
objects, many of which are known to be merging systems. 
Moreover, the dust models contain a large number of free 
parameters, which must be tuned by fitting a chosen set of 
observations, and may or may not be constant from galaxy 
to galaxy or over cosmic time. 

An alternative approach is to develop an ana- 
lytic or semi-analytic model to estimate the fraction of 
starlight that is absorbed by dust in a given galaxy, 
based on its geome try, metal content and stellar popula- 
tions. Early SAMs jCuiderdoni fc Rocca-Volmerangelll987l : 



^ The convergence is not as good in the NIR, where there are 
still significant uncertainties regarding the importance of con- 
tributions from Thermally Pulsa ting Asymptotic Giant Branch 
(TP-AGB) stars (jMar^t^Hooi). 
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Lacev et al.lll993l: iGuiderdoni et al] | l998l: iKauffmann et al 



19981 : ISomerville fc Primacklll999l : iDevriendt fc Guiderdoni 



200(3 ') approached this by assuming that the face-on B or 



V-band optical depth of the dust in the disc is proportional 
to the column density of metals in the gas phase, that the 
inclination dependence is that predicted by a simple "slab" 
model in which the dust and stars are uniformly mixed, and 
that the wavelength dependence is given by a fixed "attenu- 
ation law" , such as a Galactic extinc t ion law or the starburst 
attenuation law of ICalzetti et aD (|200(]| ). ICharlot fc Faill 
|2ooa i proposed a two-component model that separately ac- 
counts for the extinction due to diffuse "cirrus" in the disc 
and that d ue to the dense "birth clou ds" surrounding newly 
born stars. Lucia fc Blaizod l|2007l ) combined the two ap- 
proaches, using a "sla b" model to treat the cirrus compo- 
nent and adopting the lCharlot fc Falll (l2000l) model to treat 
the young stars ( ^ 10^ yr). Fontanot et al.l ()2009bl ) tested 
a wide range of such simple analytic approaches from the 
literature agai nst full radiative tra nsfer, applied within the 
MORGANA (jMonaco et alj l2007l ) SAM. They concluded 
that bulk properties (such as UV, optical, and NIR lumi- 
nosity functions) predicted by the SAMs using analytic dust 
recipes agreed quite well with the results of the full radiative 
transfer, at a fraction of the computational cost. 

With an estimate for the total energy absorbed by 
dust in hand, and assuming that all of this energy is 
re-radiated in the IR, one can use observationally de- 
rived or observationally calibrated template SEDs de- 
scribing the wavelength dependence of the dust emission 
JPeyriendt et al.lll999l: ICharv fc Elbajl200ll: JPale fc Heloul 



l2002l : iLagache et all l2004l: iRieke et al.ll2009l ). or modified 

Planck functions feaviani et al. 2003) to compute IR lu- 

1 ' 1) m " — 1 

minositi es (IGuider doni ct al. 1998; Devrie ndt fc Guiderdonil 

[200a : .Hatton et al...2003 : Blaizot et al. ,2003'). Observation- 
ally, it is known that the mid- to far-IR colours (i.e. the 
ratio of warm to cool dust) are correlated with the total 
IR luminosity of the galaxy (|Sanders fc Mirabellll99d ). Ac- 
cordingly, models based on this approach use an SED library 
indexed by the total IR luminosity; i.e., the total IR lumi- 
nosity of the model galaxy is used to s elect the appropriate 
FIR template. 'Fontanot fc Som ervillj (|20ld ) compared this 
kind of approach, again applied to the MORGANA SAM, 
with the results of coupling the s ame SAM with the full RT 
model GRASIL (|Silva et al.lll998l '). Again, the agreement for 
statistical quantities such as luminosity functions was quite 
good. 

Our goal in this paper is to develop fully semi-analytic 
models of galaxy formation that can predict photometric 
properties of galaxies from the far-UV to the far-IR with rea- 
son able accuracy. To do t his, we adopt a modified version of 
the lCharlot fc Falll Hooo') model to estimate how much light 
is absorbed by dust in each galaxy, assume that all of this 
energy is re-radiated by dust, and employ observationally- 
calibrated templates to estimate luminosities in the mid- to 
far-IR. We first ask how successfully this simple approach 
can reproduce a compilation of observations spanning the 
FIR to the FUV and a redshift range of zero to five. We 
expect such a naiVe approach to fail in some respects, and 
we attempt to identify the physical lessons that we can draw 
from the points of failure. To aid in this, we vary some of the 
more uncertain ingredients of our models to study the effect 
on the observables. The resulting fiducial models will be used 



to create m ock catalogs for pan-chromatic surveys su ch as 
CANDELS (|Grogin et al.l201ll : lKoekemoer et al.r201ll ). and 

to make predictions of the Extragalactic Background Light 
(EBL) and its build-up over cosmic time. The implications 
for gamma ray observations will be explored in a compan- 
ion paper (Gilmore, Somerville, Primack & Dominguez 2011; 
hereafter GSPD). 

This paper is structured as follows. In ij2]we describe the 
ingredients of the semi-analytic model, including our treat- 
ment of dust attenuation and emission. In [|3]we present pre- 
dictions for luminosity functions, counts, and related quanti- 
ties from the FUV to the sub-mm for several model variants, 
and compare them with an extensive compilation of obser- 
vations. We discuss and summarize our results in ^ 



2 MODELS 

2.1 Semi-analytic Models 

2.1.1 Galaxy Formation 

The semi-a nalytic models used h e re ha v e been described 
in det ail in ISomerville fc Primaclj (Il999l l. ISomerville et aP 
(|200ll) and especially ISomerville et al. (|2008al . hereafter 
SOS), and we refer the reader to those papers for details. 
Here we provide a brief summary of the basic ingredients of 
the SAM, which include the growth of structure in the dark 
matter component in a hierarchical clustering framework, 
radiative cooling of gas, star formation, supernova feedback, 
AGN feedback, galaxy merging within dark matter haloes, 
metal enrichment of the interstellar medium (ISM) and in- 
tracluster medium (ICM), and the evolution of stellar pop- 
ulations. 

We assume a sta ndard ACDM universe and a Chabrier 
IMF (|Chabrieill2003 ). We consider two sets of cosmological 
parameters here. One is the "concordance" cosmology (C- 
ACDM), Q.m = 0.3, = 0.7, Ho = 70.0, and ag = 0.9, 
which was used by SOS and has also been used by a large 
number of other studies in the literature. We also present 
results for an updated set of cosmological parameters that 
are consistent with the five year WMAP results (WMAP5): 

= 0.2S, = 0.72 , Ho = 70.0, as = O.Sl, and = 0.96 
(|Komatsu et al.ll2009l ). We note that these values are gener- 
ally consistent with those obtain ed from the analysis of the 
seven-year WMAP data release (|Komatsu et al. II2OIO'). The 
adopted baryon fraction is 0.1658. We consider WMAP5 to 
be our "fiducial" model and present results from C-ACDM 
for ease of comparison with previous work. 

The merging histories (or merger trees) of dark 
matter haloes are constructed based on the Extended 
Press-Schechter formalism using the method described in 
ISomerville fc KolattI (|l999l ). with improvements described 
in SOS. These merger trees record the growth of dark mat- 
ter haloes via merging and accretion, with each "branch" 
representing a merger of two or more haloes. We construct 
grids of "root halos" spanning the range Vvir = 30 — 1200 
km/s, where Vvir is the circular velocity at the virial ra- 
dius, and resolve the merger history of each root halo down 
to progenitors at least 0.01 times the mass of the root. For 
root halos more massive than 10^^ M0, we follow the merger 
histories down to a minimum progenitor mass of 10^" Mq . 
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We have checked that our results are robust to the chosen 
mass resolution of the trees. 

Whenever dark matter haloes merge, the central galaxy 
of the largest progenitor becomes the new central galaxy, 
and all others become 'satellites'. Satellite galaxies lose an- 
gular momentum due to dynamical friction as they orbit and 
may eventually merge with the central galaxy. To estimate 
this merger ti mescale we use a variant of th e Chandrasekhar 
formula from iBovlan-Kolchin et al.l (|2008l ). Tidal stripping 
and destruction of satellites are also included as described in 
SOS. We have checked that the resulting mass function and 
radial distribution of satellites (sub-haloes) agrees with the 
results of high-resolu tion N-body simula tions that explicitly 
follow sub-structure (|Macci6 et al.|[201oh . 

Before the Universe is reionised, each halo contains 
a mass of hot gas equal to the universal baryon fraction 
times the virial mass of the halo. After reionisation, the 
photo-ionising background suppresses the collapse of gas 
into low-mass haloes. We use the resuhs of lGnedinI \2QQ(h 
and iKravtsov et al.l (|2004l ') to model the fraction of baryons 
that can collapse into haloes of a given mass after reion- 
isation, assuming that the universe was fully reionized by 
z = 11. 

When a dark matter halo collapses, or experiences a 
merger that at least doubles the mass of the largest progen- 
itor, the hot gas is shock-heated to the virial temperature 
of the new halo. This radiating gas then gradually cools 
and collapses. The cooling rate is estimated using a simple 
spherically symmetric model, based on the following picture. 
Assuming that the density profile of the gas decreases mono- 
tonically with increasing radius, and the cooling rate is more 
rapid for dense gas, at any moment we can define the "cool- 
ing radius" as the radius within which all the gas will have 
had time to cool within a time tcooi- Then, assuming that 
the initial density profile of the gas is a singular isothermal 
sphere (pgas oc r~^), the cooling rate is given by: 

1 ^cool 1 / ^ \ 

jricooi — -^mhot 7 , (IJ 

^ ^vir ^cool 

where mhot is the mass of the hot halo g the virial 

radius of the dark matter halo, and TcooI is the cooling ra- 
dius. We calculate the cooling radius using the metallicity 
dependent atomic cooling curves of [Sutherland fc Dopital 
j 19931 ). Previous studies have used different values for the 
cooling time tcooi (e.g., the Hubble time, the time since the 
last halo merger, or the halo dynamical time). In the mod- 
els of SOS, and also here, it is assumed to be equal to the 
halo dynamical time, tdyn oc rvir/K-ir, where Vvir is the virial 
velocity of the halo. 

In some cases the cooling radius can be formally larger 
than the virial radius. In this case, the cooling rate is limited 
by the infall rate: 

T^cooi = ^nihot—^- (2) 

The cooling radius limited regime (rcooi < ''vir) is often as- 
sociated with what have been termed "hot fiows" in hy- 
drodynamic simulations, in which gas is shock heated in a 
diffuse halo and then cools. The infall limited cooling regime 
(^cooi > Tvir) is associated with "cold flows", in which gas 
streams into th e halo alon g dense filaments, without ever 
getting heated jBirnboim fc Dekelll2003l : lOekd fc BirnboimI 
l2006l : lKeres et al.ll2005l ). 



In the present SAM, we assume that the cold gas is 
accreted only by the central galaxy of the halo, but in re- 
ality satellite galaxies should also receive some measure of 
new cold gas. In addition, we assume that all newly cooling 
gas initially collapses to form a rotationally supported disc. 
The scale radius of the disc is computed based on the initial 
angular momentum of the gas and the halo profile, assum- 
ing that angular momentum is conserved and that the self- 
gravity of the collapsing baryons c auses contraction of the 
matter in t he inner part of th e halo (Bl umenthal et aDilQSd : 
iFlores et al . 1993; Mo et al.l ll998). This approach has been 
shown to reproduce the obser ved size versus stellar m ass 
relation for discs from z ~ 0-2 (|Somerville et al.ll200Sbl) . 

Star formation occurs in two modes, a "quiescent" mode 
in isolated discs, and a merger-driven "starburst" mode. 
Star formation in isolated discs i s modelled using the empir- 
ical Schmidt-Kennicutt relation (|Kennicuttlll99Sl ). assummg 
that only gas above a fixed critical surface density is eligi- 
ble to form stars. The efficiency and timescale of the merger 
driven burst mode is a function of the merger mass ratio and 
the gas fractions of the progenitors, and is based on the re- 
sults of hydrodynam ic simulations (|Robertson et al]|2006al : 
iHopkins et"ai]|2009al '). 

Some of the energy from supernovae and massive stars is 
assumed to be deposited in the ISM, resulting in the driving 
of a large-scale outfiow of cold gas from the galaxy. The mass 
outfiow rate is proportional to the star formation rate and 
inversely proportional to the galaxy circular velocity (escape 
velocity) to the power of Orh, where cvrh — 2, as expected 
for "energy driven" winds. Some fraction of this ejected gas 
escapes from the potential of the dark matter halo, while 
some is deposited in the hot gas reservoir within the halo, 
where it becomes eligible to cool again. The fraction of gas 
that is ejected from the disc but retained in the halo versus 
ejected from the disc and halo is a function of the halo circu- 
lar velocity (see SOS for details), such that low-mass haloes 
lose a larger fraction of their gas. 

Each generation of stars also produces heavy elements, 
and chemical enrichment is modelled in a simple manner 
using the instantaneous recycling approximation. For each 
parcel of new stars dm*, we also create a mass of met- 
als dMz = ydm,, which we assume to be instantaneously 
mixed with the cold gas in the disc. The yield y is assumed 
to be constant, and is treated as a free parameter. When gas 
is removed from the disc by supernova driven winds as de- 
scribed above, a corresponding proportion of metals is also 
removed and deposited either in the hot gas or outside the 
halo, following the same proportions as the ejected gas. 

Mergers are assumed to remove angular momentum 
from the disc stars and to build up a spheriod. The effi- 
ciency of disc destruction and spheroid growth is a function 
of progenitor gas fraction and merger mass ratio, and is pa- 
rameteri zed based on hydrodyn amic simulations of disc-disc 
mergers (|Hopk ins et al]|2009al ). These simulations indicate 
that more "major" (closer to equal mass ratio) and more 
gas-poor mergers are more efficient at removing angular mo- 
mentum, destroying discs, and building spheroids. Note that 
the treatment of spheroid formation in merger s used here has 
been updated relative to SOS as described in lHopkins et al] 
(2009b) . The updated model produces good agreement with 
the observed fraction of disc vs. spheroid dominated galaxies 
as a function of stellar mass. 
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In addition, mergers drive gas into galactic nuclei, fu- 
eling black hole growth. Every galaxy is born with a small 
"seed" black hole (typically ~ 100 M0 in our standard mod- 
els). Following a merger, any pre-existing black holes are as- 
sumed to merge fairly quickly, and the resulting hole grows 
at its Eddington rate until the energy being deposited into 
the ISM in the central region of the galaxy is sufHcient to sig- 
nificantly offset and eventually halt accretion via a pressure- 
driven outflow. This results in self-regulated accretion that 
leaves behind black holes that naturally obey the observed 
correlation between BH mass and s pheroid mass or velocity 
dispersion (iDi Matteo et all l2005l : [Robertson et al.l l2006bl : 
ISomerville et al.ll2008al v 

There is a second mode of black hole growth, termed 
"radio mode", that is thought to be associated with pow- 
erful jets observed at radio frequencies. In contrast to the 
merger-triggered mode of BH growth described above (some- 
times called "bright mode" or "quasar mode"), in which the 
BH accretion is fueled by cold gas in the nucleus, here, hot 
halo gas is assumed t o be accreted according to the Bondi- 
Hoyle approximation l|Bondi|[l952l ). This leads to accretion 
rates that are typically only about < 10~^ times the Ed- 
dington rate, so that most of the BH's mass is acquired dur- 
ing episodes of "bright mode" accretion. However, the radio 
jets are assumed to couple very efficiently with the hot halo 
gas, and to provide a heating term that can partially or 
completely offset cooling during the "hot flow" mode (we 
assume that the jets cannot couple efficiently to the cold, 
dense gas in the infall-limited or cold flow regime). This 
"radio mode feedback" appears to be able to quite success- 
fully solve many of the problems experienced by previous 
generations of ACDM-based galaxy formation models; for 
example, the excess number density and overly high spe- 
cific sta^JonnMion^a^esaBd^ of massive galax- 
ies (ICroton et al. l l2006l : lBowe~ et al. I I2OO6I : ISomerville et al.l 
bOOSal ). 

For each galaxy, we store a two-dimensional grid of the 
mass of stars with a given age and metallicity, separately for 
the disc and spheroid. We then convolve this distribution 
with the predictions of the simple st ellar population (SSP) 
models of lBruzual fc CharlotI (|2003l ') to obtain the SED of 
the unattenuated starlight. We use the models based on the 
Padova 1994 isochrones with a lChabrieJ l|2003l ) IMF. 



2.1.2 Dust Attenuation 

Our mod el for dus t extinction is based on the model pro- 
posed bv lCharlot fc Fall (,2000, '). As in their model, we con- 
sider extinction by two components, one due to the dif- 
fuse dust in the disc and another associated with the dense 
'birth clouds' surrounding young star forming regions. The 
1^-band, face-on extinction optical depth of the diffuse dust 
is given by 



model; i.e. the extinction in the l/-band for a galaxy with 
inclination i is given by: 



TVfi ~ Tdust.O Zcold '71cold/(''gas) 



(3) 



where rdust,o is a free parameter, Zcoid is the metallicity 
of the cold gas, mcoW is the mass of the cold gas in the 
disc, and rgas is the radius of the cold gas disc, which is 
assumed to be a fixed multiple of the stellar scale length (see 
SOS). To compute the actual extinction we assign a random 
inclination to each disc galaxy and use a standard 'slab' 



-2.5 log 1 



- exp[~Tv.o/ cos(i)] 



Tv,o/ cos(i) 



(4) 



Additionally, stars younger than tac are enshrouded 
in a cloud of dust with optical depth tbc.v = Mbctv.o, 
where we treat tuc and fiBC as free parameters. Finally, 
to extend the extinction estimate to other wavebands, we 
assume a starburst attenuation curve tea lzetti et~al]|2000l ) 
for the diffuse dust component and a power-law extinction 
curve PC (A/5500A )", with n = 0.7, for the birth clouds 
(|Charlot fc Fallll2000l ). 



2.1.3 Dust Emission 

Using the approach described above, we can compute the to- 
tal fraction of the energy emitted by stars that is absorbed 
by dust, over all wavelengths, for each galaxy. We then as- 
sume that all of this absorbed energy is re-radiated in the 
infra-red (we neglect scattering), and thereby compute the 
total IR luminosity of the galaxy Lib.. We make use of dust 
emission templates to determine the SED of the dust emis- 
sion, based on the hypothesis that the shape of the dust SED 
is well-correlated with Lir. The underlying physical notion 
is that the distribution of dust temperatures is set by the 
intensity of the local radiation field; thus more luminous or 
actively star forming galaxies should h ave a larger propor- 
tion o f warm dust, as is in fact observed l|Sanders fc Mirabell 
Il99d) . 

There are two basic kinds of approaches for construct- 
ing these sorts of templates. The first is to use a dust 
model along with either numerical or analytic solutions 
to the standard RT equations to create a library of tem- 
plates, calibrated by compari son with local protot ypes. This 
approach was pioneered by 'Des ert et al.l (Il990h. a nd has 
been followed by many other w orkers (e.g. G uiderdoni et al.l 
19981: iDevriendt et alll999l : lOevriendt fc Guiderdonill2000l ). 
Desert et al.l (|l990h posited three main sources of dust emis- 



sion: polycyclic aromatic hydrocarbons (PAHs), very small 
grains and big grains. The grains are composed of graphite 
and silicates, with small and big grains probably dominated 
by graphite and silicate respectively. The thermal properties 
of each species are determined by the size distribution and 
thermal state. Big grains are assumed to be in near ther- 
mal equilibrium, and their emission can be modeled as a 
modified black-body spectrum. However, small grains and 
PAHs are probably in a state that is intermediate between 
thermal equilibrium and single photon heating. They are 
therefore subject to temperature fluctuations and their emis- 
sion spectra are much broader than a modifled black-body 
spectrum. The detailed size distributions are modeled us- 
ing free parameters, which are calibrated by requiring the 
model to flt a set of observational constraints, such as the 
extinction/attenuation curves, observed IR colours and the 
IR spectra of loc al galaxies. Here we ma ke use of the tem- 

hereafter DGS99) 
(|l990l ). 



plates derived bv IDevriendt et al.l (ll99E 



using a similar approach to lDesert et al. 



The second approach is to make direct use of observed 
SED s for a set of p rototy pe galaxies (e.g . Charv fc Elbaj 
I2OOII : iDale fc Heloui 12003 : iLagache et all 12004 ). We also 
make use of the empirical SED templates recently published 
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Figure 1. Comparison of the dust emission templates of iRieke et alj ||2009|) (red) and iDevriendt et alj lll999l) (blue). The four panels 
show the dust emission templates used in this work for bolometric IR luminosities of 10^'^ Lq, 10^^ Lq (a LIRG), lO^'^ Lq (a ULIRG), 
and 10^^ Lq (an extremely IR-bright 'Hyper— LIRG'). 



bv lRieke et~al] (|2009l . hereafter R09). They constructed de- 
tailed SEDs from published ISO, IRAS and NICMOS data 
as well as previously unpublished IRAC, MIPS and IRS 
observations. They modeled the far infrared SEDs assum- 
ing a single blackbody with wavelength-dependent emissiv- 
ity. The R09 library includes fourteen SEDs covering the 
5.6x10^1/0 < LiR < 1O"L0 range. Examples of the DGS99 
and R09 SED templates are shown in Figure [1] 

The R09 templates have less emission in the PAH and 
mid-IR regions than those of DGS99, particularly at the 
brightest luminosities. The R09 templates are also consid- 
erably more detailed in their representation of PAH emis- 
sion. Being observationally based, the shortest wavelengths 



predicted by the R09 templates are contaminated by emis- 
sion from direct starlight. We have attempted to remove 
this component by subtracting from each template the av- 
erage amount of starlight in the SED for galaxies of that 
IR luminosity in the local universe. The R09 templates also 
end abruptly at 5 fim, and we have smoothed the transition 
to the shorter wavelength starlight regime by extrapolating 
using a power-law of slope ~ A^. Our results are not very 
sensitive to the choice of this power-law slope. 



© 0000 RAS, MNRAS 000, 000-000 



Galaxy Properties from the UV to the FIR 7 




Figure 2. Left: Galaxy luminosity in the IR relative to the luminosity in the UV vs. bolo metric luminosity. Open symbols with error bars 
are observatio nal e s timat es of this relationship for nearby galaxies from lBuat et al] ||2007|). The long dashed green line is the observational 
relation from iBelll l|2003l '). and the dotted line is the observational relation from IXu et al ] l|2006l ). Right: Galaxy lumi nosity in the IR 
relativ e to the UV vs. metallicity of the cold gas in the galaxy. The da shed (green) line is the observational relation of lHeckman et ahl 
l|l998l l and the solid green line is the relation of ICortese et al] ll2006l ). In both panels, solid black squares show the medians for our 
fiducial model, and dashed black lines show the 16th and 84th percentiles. 



2.1.4 Galaxy Formation Parameters 

The galaxy formation models contain a number of free pa- 
rameters which are tuned using observational constraints. A 
full list of the physical parameters in the semi-analytic model 
is given in SOS (Table 2). The most important parameters 
for the quantities presented in this paper are those control- 
ling supernova and AGN feedback; specifically, the efficiency 
of supernova-driven outflows (esn) ^nd its dependence on 
galaxy circular velocity (orh), and the efficiency of heating 
by "radio mode" AGN feedback (Kradio). As in SOS, the val- 
ues of these parameters are adjusted in order to obtain a 
good match to the observed z ~ stellar mass function. 
The low-mass end of the stellar mass function is insensitive 
to the AGN feedback recipe and is mainly controlled by the 
supernova feedback recipe, while the reverse is true for the 
high-mass end. The effective yield used for chemical evolu- 
tion is fixed by matching the zero-point of the galaxy stellar 
mass- metallicity relation (see SOS). The values of the param- 
eters used in the C-ACDM model presented in this paper are 
identical to those for the C-ACDM model in SOS, except that 
we have slightly increased the strength of the supernovae 
feedback (eg^ = 1.5 and arh = 2.5), and decreased the ef- 
ficiency of the radio mode feedback (/tradio = 2.5 x 10~^.) 
For the fiducial WMAP5 model presented here, we have ad- 
justed the galaxy formation parameters slightly to account 
for the modified cosmology, but they are quite similar to the 
parameters for the WMAP3 model presented in SOS. 

2.1.5 Dust Parameters 

We also have three additional parameters that control the 
dust attenuation in our model: the normalization of the face- 
on V-band optical depth Tdust.o, the opacity of the birth- 



clouds relative to the cirrus component fiBC, and the time 
that newly born stars spend enshrouded in their birthclouds, 
tsc. We first set rdust.o by matching the normalization of the 
observed relationship between I/dust /iuv vs. bolometric lu- 
minosity Lboi, where I/dust is the total luminosity absorbed 
by dust and re-emitted in the mid- to far-IR and Luv is the 
luminosity in the Far UV (~ 1500^4). The predicted median 
relation and 1-a scatter is shown in Figure [21 along with 
several observational estimates. We also show the predicted 
relationship between the cold gas metallicity and Ldust/^-uv 
compared with observations in Figured) We find good agree- 
ment with Tdust.o = 0.2. With this value, we also obtain good 
agreement with the observed optical through NIR luminos- 
ity functions at 2: = (see Section [311]). 

The birthcloud parameters /iBC and tec mainly con- 
trol the attenuation of UV light relative to longer wave- 
lengths. At 2 = 0, the g through is'-band luminosity func- 
tions are insensitive to the birthcloud parameters, while the 
FUV through u-bands are quite dependent on them. We ad- 
just these parameters in order to match the z = FUV and 
NUV observed luminosity functions, finding good agreement 
with fiBC = 4.9 and tsc = 2 x 10^ yr. 

It remains an open question whether the properties of 
interstellar dust have evolved over cosmic time and, if so, 
what impact this might have on the appearance of high red- 
shift galaxies. Recent work has suggested that, at fixed Lboi, 
galaxies may be less extinguished at high redshift than one 
would expect if the relat ion shown in Figti r e [2] w ere con- 
stant over all times (e.g. iReddv et all I2OO6I . [201 Oh . In ad- 
dition, we find, in agreement with some previous studies 
(|Lo Faro et all 12009 : Guo & White 2009;), that if we keep 
the dust parameters fixed at constant values, we under- 
predict the number of UV-bright galaxies at high redshift. 
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Therefore we introduce an ad hoc redshift dependence into 
the dust parameters, which we use in our fiducial model and 
most of its variants: rdust,o(-z) = Tdusst.oC^ ~ 0)/(l + 
both /iBC and tsc scale with above z — 1. We chose 
this dependence because it allows us to fit the bright end 
UV and B-band luminosity functions at all redshifts where 
they are well constrained observationally. We also show re- 
sults from a model with dust parameters that do not vary 
with redshift ("fixed dust" model). 



2.1.6 Model Variants 

We vary several of the uncertain ingredients of our models 
in order to study the sensitivity of our results to these as- 
sumptions. As we have already discussed, we consider two 
cosmologies, the "concordance" (C-ACDM) or WMAPl cos- 
mology and the currently favored WMAP5 cosmology. We 
consider a model variant in which instead of using the com- 
posite du st attenuation rnodel ( with cirrus plus birth clouds) 
based on ICharlot fc Falil d 2000), we instead use a fixed at- 



tenuation curve from Calzetti et al. 1 120001). We consider two 
different sets of empirical dust emission SED templates (R09 
and DGS99). As discussed above, we also consider a model 
in which the dust parameters do not evolve with redshift. 
This results in five models, which are summarized in Ta- 
ble [1] This table also specifies the line style that will be 
used for each model in the plots of the following section. 
For some quantities, some of the models produce the same 
predictions, in which case we show a subset of the models. 



3 RESULTS 

3.1 Luminosity Functions and Luminosity 
Densities 

Although our models are very similar to those presented in 
SOS, the resolution used here is much higher than the resolu- 
tion used in the simulations presented in SOS. In SOS, halos 
were followed down to a minimum mass of 10^° Mq, while 
here, root halos at z = are followed down to a minimum 
progenitor mass of 1.3 x 10* M©. Therefore we first inves- 
tigate whether changing the resolution of the merger trees 
has any effect on our results. In Fig. [3] we show the two main 
quantities that we use to normalize our models, the stellar 
mass function and the gas fraction of disk-dominated galax- 
ies as a function of stellar mass, at 2: = 0. In SOS, we stated 
that galaxies with stellar masses above ~ 10^ M© should be 
reliably resolved. We see from Fig. that the results above 
this mass change negligibly when we increase the mass reso- 
lution by almost two orders of magnitude. In addition, with 
no tuning, our model fits the observed stellar mass func- 
tion extremely well down to the smallest masses where it 
has been measured. The gas fractions also appear to match 
observations down to very low mass objects (~ lO^M©). 

Next, as a further orientation to our models, we show 
in Fig. [4] the global star formation rate density and stellar 
mass density across cosmological time. The global star for- 
mation histories are similar in the two models below redshift 
two, but the C-ACDM model has more early star formation 
because of the larger amount of small-scale power. We can 
see that there is a very small difference between the fiducial 



model at the resolution of SOS and with the higher reso- 
lution. This shows that the resolution adopted by SOS was 
sufficient to resolve the galaxies that contribute the bulk of 
the global star formation rate and stellar mass density from 
2; ~ 6 to the present. 

As a first test we investigate the luminosity functions 
at z = from the rest-UV to the IR. In Figure [5] we show 
the FUV and NUV luminosity functions from GALEX, ugriz 
LFs from SDSS, the K-band LF from 2MASS, and the total 
IR luminosity function (references given in figure caption). 
Our fiducial model agrees very well with all of these data, 
other than slightly overpredicting the numbers of the faintest 
galaxies in the GALEX FUV and NUV bands. We see that 
the model with a fixed attenuation curve (Calzetti model) is 
unable to simultaneously reproduce the UV through optical 
bands, while the composite (modified Chariot & Fall model) 
does this well. In this model, the young stars that produce 
most of the UV light are heavily enshrouded in dense birth 
clouds, resulting in an effectively age-dependent extinction 
curve. 

Next we explore the evolution over cosmic time of lumi- 
nosity functions in several wavebands. All luminosity func- 
tions are shown in the rest-frame at the indicated red- 
shift. We assume that the observations have been properly 
completeness corrected and make no attempt to determine 
whether our model galaxies would in fact obey the secondary 
selection criteria of any given survey (for example, the colour 
criteria in Lyman- break selected samples). We defer these 
sorts of more detailed comparisons to future work. Figure [B] 
shows the rest-frame ~ 1500 A LFs from z — 0.5 to z = 5. 
We see that the models produce more than enough UV lumi- 
nous galaxies at all redshifts, requiring some extinction even 
at z ~ 5. However, the model with fixed dust parameters 
(normalized a,t z — 0) systematically underproduces UV- 
luminous galaxies by a larger and larger factor as redshift 
increases. This has been found before by other studies us- 
ing semi-analytic models (|Lo Faro et al.ll2009l : lGuo fc White! 
2009). Moreover, as already mentioned, there is direct obser- 
vational evidence that high redshift galaxies m ay be less ex- 
tingu ished than their low redshift counterparts l|Reddv et al.l 
I2OIOI ). With the simple evolving dust model, we obtain fairly 
good agreement over all redshifts, although we somewhat 
overproduce low- luminosity galaxies at high redshift. Obvi- 
ously, we could have adopted a more complex dust model 
with luminosity dependent evolution in the dust parameters 
to get a better overall fit. Alternatively, this could be an in- 
dication tlM^Jow^^nMssealajd^ are forming too early in the 
models (|Fontanot et al. I l2009al ). 

Figure [7] shows a similar comparison in the rest B-band. 
We see a similar discrepancy with our non-evolving dust 
model to that observed in the UV: we need to adopt lower 
extinctions in high redshift galaxies to reproduce the num- 
ber of luminous galaxies. The discrepancy is smaller, and 
sets in at higher redshift, but is still pronounced by z ~ 3. 
The combination of the rest UV and optical observations are 
what forced us to adopt evolving parameters in both com- 
ponents of our two-component model (the cirrus and the 
birthclouds). 

Next, in Figure [H] we show the comparison with the 
rest-frame K-band luminosity functions in the same redshift 
bins. Here, interestingly, the models over-predict the num- 
ber of galaxies at high redshift, particularly below the knee 
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Table 1. Summary of Models 



cosmology dust attenuation 



dust emission 



dust parameters designated line style 



WMAP5 composite Rieke et al. (2009) evolving 

WMAP5 composite Rieke et al. (2009) fixed 

WMAP5 Calzetti Rieke et al. (2009) fixed 

WMAP5 composite Devriendt et al. (1999) evolving 

C-ACDM composite Rieke et al. (2009) evolving 



WMAP5 fiducial 
WMAP5+fixed dust 
WMAP5+Calzetti 
WMAP5+DGS 
C-ACDM 



solid black 
dash-dotted purple 
dashed red 
long-dashed blue 
dotted black 




Figure 3. Left: Stellar mass function at ^ = 0. The solid line shows the fiducial (WMAP5) model, dotted line shows the C-ACDM model, 
and the dashed gray line shows the fiducial mode l at th e resolution of the simulations pres ented in SOS. Sy mbols show the obs ervationally 
derived stellar mass functions fro m lBaldrv et al] ||2008| . blue squares). iBaldrv et al.l ll201ll . purple crosses). iPanter et al.l | |2007| . dark green 
triangles), and lLi fc White! | |2009| . li ght green hexagons). Note that the quoted values for the observed mass functions at rristar IO^'^Mq 
are likely to be lower limits, due to surface brightness selection effects. Right: Gas fraction for central disk-dominated galaxies in the 
fiducial model (gra y shading as i n SO S). Solid and dashed lines show model median and 16 and 84th percentiles. Large open circles show 
the observations of iKannappanI ||20041. Fille d squares show g as fractions from galaxies in the THINGS survey llLerov et al.ll2008l) . and 
small filled circles show observations from iBaldrv et al.l l|2008l ) . 



in the l uminosity funct i on. Th is is consistent with the find- 
ings of iFontanot et al] (|2009al ). who showed that three in- 
dependently developed semi-analytic mode ls all overproduce 
low-m ass galaxies at high redshift (see also lMarchesini et al.l 
I2009D . The good agreement of our models at the bright- 
est K-band magnitudes also indica tes that, contrary to the 
findings of iHenrigues et ahl (|2010l ). there is no need to in- 
voke an enhanced contribution to the Near-IR due to the 
Thermally Pulsating Asymptotic Giant Branch (TP-AGB) 
stellar phase, a s in e.g. the stellar population models of 
iMarastonI (j2005h . 

We now move from bands that are dominated by hght 
directly emitted by stars to light that has been reprocessed 
by dust. Figure [9] shows the LP in the Spitzer IRAC rest 8 
/xm band. Our models produce about the right number of 
galaxies around the knee in the LP, but underproduce lumi- 
nous galaxies, especially at high redshift {z ~ 1-2). As one 
can see from Pigure [T] this portion of the spectrum is very 
complex, and highly dependent on whether the spectrum is 
dominated by emission from PAHs, or absorption features 
such as the deep Silicate absorption seen at ~ 9pim. In ad- 



dition, this part of the spectrum is particularly subject to 
contamination by obscured AGN. If there are a significant 
number of buried AGN at « ~ 2, as has been suggested 
iDaddi et al.ll2007l ). this could at least partially account for 
the discrepancy between the observations and our model. 
However, it would also be unsurprising if our simple unevolv- 
ing template approach were simply inadequate to accurately 
model the very complex physics that determines the 8 /im 
flux in galaxies. Indeed, we see a large difference between 
the predictions using the R09 templates and those using the 
DGS99 dust emission templates. Very similar remarks apply 
to the comparison between the Spitzer MIPS 24 fim LPs in 
Pigure 1101 although here the models agree reasonably well 
with the observations at z = 0.5 and 2 = 1, and show a 
somewhat milder discrepancy at z = 2. The differences in 
the predictions with the two different dust emission tem- 
plates are also considerably smaller than at 8 /im. 

Pigure [TT] shows the LP at rest 250 /xm, compared with 
early results from Herschel. At z = 0.5, the only redshift 
where the LP at this wavelength has been measured to date, 
the models overpredict the number of luminous galaxies. 
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Figure 4. Left: Global star formation history. Right; global stellar mass assembly history. Both panels: solid black lines show the 
predictions of the WMAP5 model; dotted lines show the C-ACDM model; dark gray dashed lines show the fiducial model at the 
resolution of the simulations presented by SOS. Symb ols show a compilation of observational estimates (references given in 808). The 
long-dashed lines are the observational estimates from lHopkins fc BeacomI ||2006| '). 



The last luminosity function we examine is that in Figure 
1121 for the total integrated IR from 8-1000 /xm, estimated 
from multi-wavelength observations. Here we see fairly good 
agreement between the models and the observational esti- 
mates, although with a small deficit of luminous galaxies, 
particularly at 2 ~ 2. However, considering that the true 
errors on the observations are probably considerably larger 
than the error bars shown, this level of agreement is encour- 
aging. It is also encouraging that the different model variants 
produce very similar results for this quantity. 

Figure [TH] shows the total emissivity from all galaxies 
as a function of wavelength in our fiducial model, at various 
redshifts. One can see a slight shift in the wavelength of the 
peak of the dust emission, as well as the ratio of warm to cold 
dust with redshift. Because we have assumed non-evolving 
dust emission templates in this work, this is entirely due 
to the changing mix of galaxies of different total IR lumi- 
nosities with time (i.e. the declining contribution of very IR 
luminous galaxies with decreasing redshift). Figure [T4l shows 
the integrated luminosity density as a function of redshift in 
the far and near-UV, the B-, z-, and K-band, and the to- 
tal IR, as a function of redshift for all the models presented 
in this work. The model predictions are compared with ob- 
servational estimates obtained by integrating observed LF's 
(see figure caption for references). Such observational esti- 
mates in the UV and IR bands are often used to estimate 
the global star formation rate density (e.g. as shown in Fig- 
ure |4| . It is interesting that while our fiducial model is in 
good agreement with the observations in the Far-UV over 
the whole redshift range, it is somewhat low compared with 
the total-IR luminosity at z ~ 1-2. As already discussed, 
this could be due to the total IR having been overestimated, 
or could indicate that the star formation rate in the models 
is too low. Also interesting to note is that the UV and IR 
luminosity densities peak at a higher redshift [z ~ 3) than 



the z- or K-band, because the latter arise from older stellar 
populations. 



3.2 Galaxy Counts 

Figures [15] and [16] show galaxy number counts from the UV 
through the sub-mm. This is an important cross-check on 
the luminosity function comparison as many high-redshift 
samples do not have spectroscopic redshifts available par- 
ticularly for the fainter objects. Our fiducial model provides 
quite good agreement with the data at all UV and optical 
wavelengths. There is a factor of 2-3 excess of faint galaxies 
in the models in the optical and NIR, which corresponds to 
the excess of faint galaxies at z ~ 1-2 seen in the luminos- 
ity function comparison. The agreement at 3.6 and 8 /^m is 
also quite good, except below ~ 30/iJy, where the models 
overpredict faint galaxies. The agreement is also good for 
MIPS 24 and 70 /xm, except for 24 ixm sources between 100 
/^Jy to 1 mjy. The models with the R09 templates underpro- 
duce galaxies in this range, while the model with the DGS99 
templates performs better. This is because of the much lower 
fluxes in the 10-20 /^m range in the R09 templates relative 
to the DGS99 templates. None of the models reproduce the 
"bump" seen in the observed 24 /im counts between ~ 0.1 
and 1 mJy. The agreement with the 70 ^ra counts is fairly 
good, but the observed counts do not go as deep. 

Moving towards longer wavelengths, the agreement be- 
comes increasingly poor. Our models underpredict galaxies 
fainter than ~ 80/iJy at 250 /im by a factor of 4-5, and un- 
derproduce sub-mm galaxies at 850 /^m by more than an or- 
der of magnitude. It is well-known that CDM-based models 
in general have difficulty reproducing the observed popula- 
tion of sub-m m ga laxies (e.g. Dcvricndt & Guidcrdoni 200_3; 
iBaugh et aklboOSl ). We discuss different possible interpre- 
tations of these results in Section [J] and we present more 
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Figure 5. Luminosity functions at 2 = in the FUV, NUV, SDSS ugriz and K bands, as well as the total IR (8-1000 fim). The solid 
black line is our fiducial WMAP5 model using the composite dust recipe (Chariot- Fall), and dashed red shows the model with a single 
component dust model (Calzetti). The dotted black line sh ows the C-ACDM m odel. The black long-das hed line shows th e predictions of 
the fiducial m odel without d ust attenuation. Data are from lWvder et al ] l|2005l 'l (blue points, (z^=0.05). [Bell et al.l ll2003l 'l (green points), 
and 'Rodighicr o^ et al.l l l201(Jl (red points, (z)=0.15). For the total IR panel, the x-axis shows the logarithmic luminosity in solar units, 
and the y-axis has units of N dex""*^ Mpc"'^; all other axes are as indicated. 



detailed predictions f or galaxies in the H erschel PACS and 
SPIRE wavebands in lNiemi et al] (|2012l ). 



3.3 The Extragalactic Background Light 

One of the major goals of this work is to predict the in- 
tegrated Extragalactic Background Light (EBL). The EEL 
consists of light emitted at all epochs, modified by red- 
shifting and dilution due to the expansion of the universe. 
Because the EBL consists of light emitted by all types of 



sources at all cosmological epochs, it forms a unique record 
of the history of photon production in the universe. A de- 
tailed measurement of the EBL flux can potentially inform 
us about the existence or nonexistence of photon sources be- 
yond those that can be resolved in galaxy surveys, and its 
SED encodes the details of the redshifts and spectral charac- 
teristics of these sources. As discussed in our companion pa- 
per (GSPD), the EBL also affects the propagation of gamma 
rays in the GeV and TeV bands. The EBL at a redshift zo 
and frequency vq in proper coordinates can be computed in 
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Figure 6. Luminosity functions in the rest frame 1500 A UV band at several nonlocal redshifts. The solid black line is our fiducial 
WMAP5 model using the evolving composite dust recipe (Chariot-Fall), the dash-dotted purple shows the composite dust model with 
fixed parameters, and dashed red is with the single component dust model ( Calzetti). The b l ack lo ng-dashed line shows th e predictions 
of the fiducial model wit hout dust attenuati on. Blue squares are data f romlArn outs ct alJ jjgooj , violet circles are from iHathi et alj 
1I2OIOI '). red stars are from lReddv et ah and green circles are from lBouwens ct al. ( 20071 )! 



our model as the following integral over emissivities from 
our predicted galaxy population: 

71 \ ^ r '^^ (l + ^o)' I ^w {K\ 
J{uo,zo) = - j^^ -^^—^e{.,z)dz, (5) 

where e{p, z) is the galaxy emissivity at redshift z and fre- 
quency u = i^o(l + 2)/(l + 20), and dl/dz is the cosmological 
line element, defined as 



dz (1 + z)Ho ya^TTlF+nl 



for a flat ACDM universe (|Peebleslll99l ). We assume here 
that the EBL photons evolve passively after leaving their 
source galaxies and are not affected by any further inter- 
actions except for cosmological redshifting. This is an ac- 
ceptable approximation for photons at energies below the 
Rydberg energy of 13.61 eV. 

We note that our estimate of the EBL does not include 
the contribution from hght radiated by AGN. However, pre- 
vious studies have shown that this contribution should be 
less than ^ 10 — 20% in the mid-IR, and smaller (a few per- 
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Figure 7. Rest frame lum i nosity functions in tfie B-band at several redshifts. Line types are the sam e as in Figure [6] Red squares are 
data from [Salimbeni et al. I l l2008h . and open blue hexes are from the study of iGiallo ngo et aL I 1I2OO5I '). The green shaded regions in the 
top two panels are the Schechter fits to DEEP survey luminosity functions llFaber e'tani2o'o7l ) . with la errors. Green circles in the two 
middle panels are from lMarchesini et al ] l|2003). 



cent) at other wavelengths (see iDonu'nguez et al.|[201ll . and 
references therein). 

Observational estimates of the EBL are obtained either 
by direct detection or by integration of galaxy counts. Di- 
rect detection is complicated by foreground emission from 
our own galaxy and reflected zodiacal hght from our sun, 
which are much brighter than the EBL across most of the 
optical and IR spectrum l|Hauser fc Dwelj|200ll ). Integrated 
galaxy counts provide a firm lower limit, but there has been 
considerable debate about how much light these estimates 
might be missing because of undetected, low surface bright- 
ness galaxies or the faint extended wings of detected galax- 



ies, which can be difficult to disentangle from the back- 
ground. Figure [17] shows the predictions of our five model 
variants along with a compilation of recent observational es- 
timates of the EBL at various wavelengths from both meth- 
ods. There is a significant gap between the direct detection 
measurements and the integrated counts, with the former 
always lying higher than the latter, indicating that there 
must still be biases or systematic sources of error in one or 
both methods. Unsurprisingly, since we have already seen 
that our models reproduce the galaxy counts at most wave- 
lengths, our model predictions all lie close to the estimates 
from integrated counts. It is noteworthy that although our 
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Figure 8. Rest fr ame luminosity f unctio ns in the K band at several re dshifts. Line types ar e the same as in Figure |6] Blue squares are 
observat ions fromlCirasuoIo etaL I ll2010h . red stars at 2 = 0.5 are from IPozzetti et al. and open green stars at redshifts 1 and 2 

are from lCaputi et ~ I ||2006|) . Note that all observations have been converted to the AB magnitude system. 



models, if anything, seem to ouer-produce faint galaxies, 
the integrated EBL predictions are nowhere nearly as high 
as the direct detection estimates. The model with the one- 
component (Calzetti) attenuation law over-produces the far- 
UV EBL and is low in the mid-IR. This is because in the 
two-component (modified Chariot-Fall) dust model, young 
stars are enshrouded in dense birth clouds with higher op- 
tical depth. The largest difference between the models with 
different dust emission templates is seen in the mid-IR, with 
the DGS99 models predicting a higher EBL in the mid-IR 
than the R09 models. The two models bracket the error- 
bars on the existing observations. Hopefully new observa- 
tions will tighten up these constraints. The peak of the EBL 



at ~ 100 — 250 //m is a bit low compared with observational 
estimates, but is within the errors. In terms of the overall 
partition of the EBL in the UV-NIR vs. NIR-FIR regimes, 
our models are in good agreement with the observational 
constraints (see GSPD). 



4 DISCUSSION 

4.1 Comparison with Previous Results 

iGuiderdoni et ahl l|l998l ) was one of the earliest studies using 
"forward evolution" , cosmologically motivated CDM models 
combined with modelling of both dust extinction and emis- 
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Figure 9. Rest frame luminosity functions in the IRAC 8 fim band at several redshifts. Long-dashed blue lines show the predictions of 
the model with the DGS99 dust emission templates. Other lines show results with the IR templates of R09: solid black is our fiducial 
WMAP5 model using our evolving composite dust recipe, dashed-dotted purple shows the model with fixed dust para meters, dashed re d 
show the single component dust (C alzetti) model, and dott ed black is the C-ACDM model. Green stars are data fromlPai et al 
Blue hexagons show the data from iRodighiero et al. I ||2010|'I . and red stars are measurements from lCaputi et al] | |2007|) . 



sion. This work was extended in iDevriendt &: Guiderdonil 
l|200(]|) and subsequent papers by the GALIC S cohabora- 
tion (jHatton et al] l2003l : iBlaizot et all |2004| '1. They cou- 
pled a semi-analytic model with analytic recipes for dust 
attenuation and theoretical templates for dust emission, 
very much in the sam e spiri t as our work here. The 
IDevriendt fc Guiderdoni ()2000l l models overpredicted the 
UV and optical counts, possibly because of their use of a sin- 
gle dust attenuation relation rather than a two-component 
model like our modified Chariot-Fall model. They also found 
that their standard model failed to reproduce enough sub- 



mm galaxies at 850 /xm, but were able to fit the sub-mm 
counts by introducing by hand a population of heavily dust 
extinguished starburst galaxies at high-redshift. 

The Durham group has coupled their s e mi-analytic 
model of galaxy formation (|Cole et al.l 1 19941 : I Cole et al] 
l200d ) with the dust models and radiative transfer code 
GRASIL, developed by ISilva et al.l (|l998l 'l. An important 
feature of the GRASIL approach is that the dust emis- 
sion SED is computed self-consistently based on the as- 
sumed properties of the dust and the predicted radiation 
field. This work is presented in an extensive series of papers 
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Figure 10. Luminosity functions in the rest-fra me MIPS 24 fim band. Line types are the sam e as in Figure|9l Blue hex agons show the 
observ ational data from lRodighiero et al.l l l201Clt). and re d stars are from lBabbedge et al.l l l2006f) . Green asterisks are from lMagnelli et ahl 
1I2OIII '). Orange squares are froni lRuiopaiaxri'etaL I 1I2OIOI '): we have interpolated these points from data presented for two different redshift 
bins. 



llGranato et all I2OO0I: iBaueh et al] l2005l: ISwinbank et all 



l2008l : lLacev et al.ll2008l . l20ld').lGranato et a.1.1 l|2000l ) showed 



that the fiducial model of lCole et al.l ||200C( ). when coupled 
with the GRASIL machinery, could reproduce local galaxy 
luminosity fu nctio ns from the UV to the FIR (12-100 fim). 
However. iBaugh et al.i (,2005. ) showed that this model failed 
to reproduce sufficient numbers of bright sub-mm galaxies 
by an order of magnitude or more. They therefore made sev- 
eral modifications to their model. They modified the star for- 
mation recipes in their model in two ways: 1) they adopted a 
quiescent star formation recipe with a constant star forma- 
tion timescale, rather than a standard Kennicutt-Schmidt 



relation, leading to larger gas reservoirs in high redshift 
galajcies; 2) they adopted an efficient starb urst mode in mi- 
nor as well as major mergers. As shown bv lSomerville et al.l 
(l200ll ). this leads to a much larger population of luminous 
starburst galaxies at high redshift. In addition, they adopted 
a top-heavy stellar Initial Mass Function (IMF) in bursts, in 
which the slope of the IMF dN/dlnm oc is a; = over 
the whole range 0.15 < m < 125 M0. In addition to produc- 
ing much more UV light per unit mass of stars formed, the 
top-heavy IMF also produces much larger amounts of met- 
als and dust. They found that all of these changes combined 
were needed to reproduce the sub-mm counts. The same 
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Figure 11. The rest-frame l uminosity functions at 250 /xm. Line types are the same as in Figure (9] and results from early Herschel 
observations l lDve et al.ll201Cl') at z = 0.5 are shown by the blue hexagons. 



model simultaneously reproduces the rest-UV and optical 
LFs at high redshift. 

iLacev et al.l (|2008l ') argue that the top-heavy IMF is 
needed in order to reproduce the evolution of the mid-IR 
LP as measured b y Spitzer, and th e counts in Spitzer bands 
as well. However, ISwinbank et al.l (|2008l ) showed that the 
K-band and IRAC 3-8 /xm luminosities of both SMGs and 
LBGs in the Lacey-Baugh models are a factor of ten too low 
(because of the dearth of low-mass stars). Moreover, the pre- 
dictions of the models do not seem to be in good agreement 
with early results from BLAST and Hers chel at 350-500 /im 
l|Lacev et al.ll2010l : IClements et al.ll2010l '). 

Comparison between our results and those of Lacey et 
al. is complicated because not only the approach for treat- 



ing dust is different, but also many of the other model in- 
gredients are significantly different. Perhaps most signifi- 
cantly, the published Lacey-Baugh models do not include 
"radio mode" feedback from AGN, which is now commonly 
adopted in semi-analytic models in order to pr event the for- 
mation of grossly over-massive galaxies (e.g. ICroton et al.l 
l2006l : iBower et al] l2006l : ISomerville et al.l l2008al ). In an at- 
tempt to solve t he overcooling problem via other means, 
I Cole et al] (|2000l ) adopted modified hot gas profiles that 
suppressed early cooling in massive haloes. This al so sup- 
press es the early for mation of massive gala xies ( Bower et alj 
l2006f ). Interestingly. iFontanot et al.1 l|2007l ) use basically the 
same approach to modeling dust (GRASIL) within a differ- 
ent semi-analytic model, and find that they are able to repro- 
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Figure 12. Luminosity functions for the total IR emission. Line types are the same an in Figure[9l however as the total IR predictions are 
insens itive to the templates used, the model with t he DGS99 templates is not s hown. Blue hexagons are the data from lRodighiero et al.l 
||2010| '). green stars are from lLe Floc'h et al. and red diamonds are from lCaputi et al.1 l^ol). 



duce the K-band and sub-mm counts with a Salpeter IMF. 
They ascribe this success to an improved cooling model. 
However, their models overproduce massive galaxies at low 
redshift (z < 1). Clearly, more work is needed in order to 
develop a model that reproduces all the available observa- 
tions. 

One should also keep in mind that the GRASIL 
dust+RT models used in the work described above as- 
sume a simplified, regular geometry (spheroid-|-disc) . Par- 
ticularly when considering populations that are likely to 
correspond to merger-driven starbursts, such as the SMGs, 
it is probably important to capture the complex physics 
and geometry of these systems. Some significant progress 



has been made recently in this regard, using hydrodynamic 
simulations of isolated and merging galaxies in combina- 
tion with the polyc hromatic Monte Carlo Radiative Trans- 
fer code SUNR ISE ('jonsson"2006':'jonsson et al."2006l. [2oTol : 
Ijonsson fc Prim ack 2010). Narayanan ct al. (2010a) studied 
the SEDs of merger simulations combined with SUNRISE 
RT calculations and identified objects wit h properties con- 
sistent with both bright and fainter SMGs. iNaravanan et akl 
(|2010bl ) used similar simulations to propose a physical model 
for the 2 ~ 2 DOGs (Dust Obscured Galaxies - optically 
faint galaxies identified at 24 /xm), which partially overlap 
with the SMG population. The SEDs from these simulations 
can be combined with predictions of merger rates from semi- 
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Figure 13. Integrated luminosity density as a function of wavelength in our WMAP5 fiducial model, shown at various slices in 
redshift. All data shown are measurements in t he local universe. Measurements are from GALEX (blue c ircle V SPSS (red st ars; 
iMontero-Dorta fc Pradall2009l) . 6dF (cva n squares: I Jones et aIll2006^ ■ and 2MASS ( green stars: ICole et al.ll200ll andlseir^ayliool). In 
the mid- and far-IR, the orange squares l|Soifer fc Neugebauedll99ll l and blue stars (Takeuchi et al.ll200ll 'l" show observations from IRAS 
and SCUBA. 



empirical models like those of lHopkins et al.l l|201(f ) , or with 
semi-analytic models, to compute statistics of the popula- 
tion, such as counts (Hayward et al. in prep). 

There is an extensive literature on "backwards evolu- 
tion" models for galaxy counts and the EBL, which we re- 
view in our companion paper GSPD, but do not discuss 
here. Recently, empirical and "semi-empirical" approaches 
have been adopted by several authors to predict the EBL. 
lYounger fc HopkinsI (|2011. ) used the observed stellar mass 
function at different redshifts, in combination with a semi- 
empirical model of galaxy evolution and template SEDs 



fromlCharv fc Elbaj (1200 it ) to predict the mid to FIR EBL. 
iDommguez et aL ( 2011 . Dll) made use of empirical tem- 
plate SEDs and observed fractions of 25 different galaxy 
types to predict the EBL and its evolution. An explicit com- 
parison with the luminosity density evolution estimated by 
their approach is shown in Figure 1141 The Dll estimates 
are anchored to the observed K-band luminosity functions, 
which our semi-analytic models reproduce fairly well, so the 
predictions are very similar at NIR wavelengths. At shorter 
wavelengths (optical and UV), the Dll approach predicts 
lower luminosity densities at high redshift than our SAMs. 
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Figure 14. Luminosity density as a function of redshift for various wavelengths, as well as for the total IR luminosity (8-1000 /im). Line 
types are the same as those discussed in previous plots; see captions for Figu res [6] and [9] and Table [T] Additionally, orange lines with 
alternate short and long dashes are the predictions of the empiric al model of iDommguez et al.l 1 I2OIII ). for comparison with our work . 
Observational data arc as follows: 1500 A: Blue squares are from lPahlen et al.l ll2007t ). red stars are from [Schiminovich et all l|2005h . 
green stars are fr om Bouwcns et al. (2007), and orange circles are froml Reddy et ahl fcoij^. The so lid purple circ l e is a local measurement 
with G ALEX bv lWvder et al.l 1)20051). 2 800 A: Blue sq uares and the purple circle are again from'Dahlcn et al.' (2007') and lWvder et al.l 
||2005| '). respectively. Red stars are from lGabasch et al.l ({2006) . B-band: Bl ue squares are from lOa hlen ct al. (2005), and are incomplete 
at higher redshifts. DEEP and COMBO-17 data from TPaber et al.l are sho wn as red stars a nd open red squ ares, respectively 

1 these are very similar and difficult to differentiate here) . Other data sho wn are from|Wolf et ah H2OO3!) reen star) andlMarchesini et al.l 
2OO7I ) (open purple hexes), z-b and: Local measurem ents are provided bv lMontero-Dor'tT'^r^ rada (200 3) (red) andlBlanton e t al. (20031*) 
(green). Blue squares are from | Gabasch et al" I ( l2006h . K-band: The loca l determination is from iKochanek et al. 1 (l200ll ). High redshift 
data are from lBarro et al (blu e squares) and lArnouts et al.l l|2007l ) (open r ed hexagons). Tota l IR Luminosi ty: observational 

estima tes of the IR emissi vity are fromlCaputi et ahl ( 200?! ) fopen blue pentagons) , iReddv et al.l ([200^ (green circles) , iRodiehiero et aTl 
\2m& i (purple stars), and lLe Floc'h et al.l (|2005l ') (red crosses). 
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Figure 15. Number counts in the GALEX UV bands and the four HST ACS bands. Line types are the same as in Figure |9] note that 
some models do not deviate significantly from the fiducial WMAP5+evoIving dust model (soUd black line) and ar e therefore not v isible. 
Note that results here h ave been rescaled to a Euclidean geometry. In the UV bands, data are from G ALEX iXu et al.ll2005l. green 
squares), STIS on H ST ({Gardner et aLlboOOl . purple asterisks), and the balloon-borne FOCA experiment l|lglesias-Paramo et al ] |2004l: 
iMilliard et"al]|l992l . red stars and open p entagons respectively). The FOCA points have been converted to the GALEX bands using 
the method described in IXu et al.l 1I2OO5I ). Blue crosses are from HST ACS/SBC observation s of mul tiple fields in G OODS-N and -S 
dVover ct al. 2011). In the ACS bands, red, blue and green squares are from the compilation bv lDolch fc Fcrgus oq l|201ll). which includes 
data from the Hubble Ultra-Deep Field. Additi onal data in orange fro m SDSS-DR6 are fro m iMontero-Dorta fc Pradal ||2009|). In the 
K-ban d, we show data fro m 6dF (orange crosses. ] Jones et al.ll2006il . from lKeenan et al.l 1I2OIC1I . open red hexagons), and from lBarro et al] 
||2009| . blue squares), and iMcCracken et al. 1 1 12010 " green pentagons). All observational data have been converted to AB magnitudes. 



These differences are d ue to the use by Dll of SWIRE tem- 
plates (|Polletta et al.l IToOTi ) from the UV to IR, while in 
our approach we model the star formation history and dust 
attenuation of each galaxy. In the Far-IR, Dll estimate a 
higher and sharper peak at z ~l-2 (again because of the use 
of different SED templates) , in better agreement with obser- 
vations, and a steeper decline at higher redshift 2^2. Note 



that the observed K-band luminosity functions that ground 
their empirical approach are available only up to z 4, and 
the results shown at higher redshifts are extrapolations. 
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Figure 16. Number counts from four Spitzer (IRAC and MIPS) infrared bands, as well as Herscliel 250 /xm and SCUBA 850 /im. Line 
types are tfie same as in Figure (9] for clarity models si milar to the fiducial m odel are not shown. Results are scaled to a Euclidean 
geometry. Solid bl ue circles in the 3.6 IRAC band are from lSanders et al. I1 I2OO7I) : all other points in the IRAC 3.6 a nd 8.0 bands are from 
iFazio et aP l|2004h . The MIPS data at 24 /im shown here are the S-COSMOS 'Extragalac tic Wide' points from ISanders et al.l l|2007l l 
(green hexes), and the Wide and Deep Legacy Survey points from [sithermin et al.l 1 I2OIOI) (blu e squares). At 70 /xm da ta shown are 
the no rmal (blue squares) and stacked (cyan squares) measurements from B ethermin et al.l ll 2010l). wh ile red stars are from lFraver et aP 
1I2OO6I '). Herschel data at 250 fira are from Clements ct al. (2010, red squares) and lci'en^^'aL moid, blue stars); the latter is from the 
spline model with FIRAS priors. We show data from the SCUBA SHADES survey Icoppin et al.ll2006h at 850 /xm in the lower-right 
panel. 



4.2 Summary and Conclusions 

We have presented predictions for the luminosity and flux 
distributions of galaxies from the far-UV to the far-IR and 
over the bulk of cosmic history [z = 0-6). Our predic- 
tions are based on semi-analytic models of galaxy formation, 
set within the hierarchical Cold Dark Matter paradigm of 
structure formation, and including modeling of gas cooling, 



star formation, stellar feedback, chemical enrichment, and 
AGN feedback. In addition, crucial to the present study is 
modeling of the attenuation and re-emission of starlight by 
dust in the interstellar medium of galaxies. We use a simple 
but physically motivated analytic approach to estimate the 
dust attenuation as a function of wavelength . In our fiducial 
models, based on the approach proposed bv lCharlot fc Falll 
( 2000. ) , young stars are enshrouded in dense "birth clouds" , 
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Figure 17. The predicted integrated Extragalactic Bacltground Light spectrum, compared with experimental constraints. Line types 
follow the convention of earher plots. Data: Upward pointing arrows indicate lower bounds from number counts; other symbols are results 
from direct detection experime nts. Note that some p oints have been shifted slightly for clarity. Lower li mits: The blue- v iolet triangles are 
results from Hubble and STIS ijGardner et al.ll200cfl. while th e purple open triangles are from GALEXJXu ^t al.*'200 5j). T he solid green 
and red triangles are from the Hubble Deep Field jMadau fc~ P ozzctti 2000) and Ultra Deep Field jDoIch & : Ferguson 201]J) respectively, 
combined with ground based-data, and the solid purple triangle is from a measurement by the Large Binocular Camera l|Grazian et al.l 
[2009). In t he near-IR J, H, a nd K bands, open violet stars are the limits from Kccnan et al. (2010.). Op en red triangles are from IRAC 
on Spitzer llFazio et al.ll20o"3) , and the purple tr iangle at 15 ;jm i s from ISOCAM (Hopwoo d et al.ll2010l) on ISO. The lower limits from 
MIPS at 24, 70, and 160 /tm on Spitzer are from lBethermin et al.l ll2010l '), ICharv et al. ( 2004 ) . iFraver et al.l | |2006|'I. andlDole et al. 1 1I2OO6I I 
(solid blue, solid gold, open gold, and open green, respectively). Lower limits from Herschel nu mber counts llBerta et al.ll2010l) are shown 
as solid red triangles. In the submillimeter, limits are presented from the BLAST experiment jDevlin et_aL 2009lV Direct D etection: 
In the optical, orange hexag ons are based on da ta from the Pioneer 10/11 Imaging Photopolarimeter l lMatsuoka et al.ll201l] '). The blue 
star is a determination from lMattila et"aL 1 ll201ll'l . and the triangle at 520 nm is an up per limi t from the same. In the ne ar-IR, the points 
at 1.25, 2.2, and 3.5 iJ, ra are based upon D IRBE data with fore ground subtra ction: IWrighj l|200ll . dark red squares'). | C arnbres y et al.l 
1I2OOII . orang e crosses), iLe vcnson & Wright | |200S| . red diamond). Icoriian et ahl 1 I2OOOI . purple open hexes) , IWright fc Reesd I2OOO'. green 
square), and [Levenson et al. (2007, red asterisks). In the far-IR, di rect detection mea surements are shown from DIRBE (jWright 20o3; 
ISchlegel et al.lll99Sl . blue stars and solid red circles), and FIRAS jFixsen et al.lll998l , purple bars). Blue- violet open squares are from 
direct photometry with the AKARI satellite l lMatsuura et al1l2010l) . 



while older stellar populations are embedded within a more 
diffuse "cirrus" component. Stars emerge from the dense 
birth clouds as they age. This two-component dust model 
results in an effectively age-dependent attenuation relation, 
such that younger stars are more extinguished. We find that 
the two-component model gives much better agreement with 
the UV-optical colours of galaxies than the widely used ap- 
proach of a fixed attenuation curve. 

We then assume that all light absorbed by dust is re- 



radiated in the IR, and use a set of "template" SEDs to 
estimate IR luminosities. Our current assumption is that 
the total IR luminosity of a galaxy determines which tem- 
plate SED to use, based on the empirical correlation between 
bolometric or total IR lum inosity and dust tempera ture in 
local LIRGS and ULIRGS (|Sanders fc Mirabelll 19961 ). How- 
ever, this is certainly too simplistic, and a goal of our future 
research is to try to understand and characterize how the 
physical parameters of galaxies impact the shape of their 
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FIR SEDs. One avenue towards this goal is to use detailed 
dust and radiative trans fer models implemented within hy- 
drodynamic simulation s (|jorissonll200 9; iJonsson et al. I l2006l . 
I2OI0I : [Narayanan et al ] |2010al lbl). Another approach is to use 
the much richer set of mid- and FIR data that is becoming 
available, spanning a broader range of cosmic epoch and 
galaxy typ e, to develop a mor e complete set of template 
SEDs (e.g. lCharv fc Popdl2010l l. 

We found that in order to fit the luminosity functions 
and counts in the UV and (to a lesser extent) optical, it 
was necessary to adopt a dust optical depth normalization 
that varied with redshift. This has the net effect that galax- 
ies of a given bolometric luminosity are less extinguished 
at high redshift, and could be interpreted as an evolving 
dust-to-metal ratio or as a geometric effect (e.g. perhaps 
the distribution of gas and dust relative to the stars is dif- 
ferent in high redshift galaxies). This result has also been 
found by othe r groups working with semi-analytic models 
ijLo Faro et al . 2009; Guo & White 2009) and is sup ported 
by direct observational evidence (jReddv et al.ll2010l ). How- 
ever, our current approach is completely ad hoc (we simply 
tuned the dust normalization parameters to match the ob- 
served FUV luminosity functions from z — 0-5) and it would 
be desirable to develop better observational constraints as 
well as a deeper physical understanding of this effect. With 
the evolving dust model, we find very good agreement with 
observed far-UV luminosity functions from z ~ 0-5 and B- 
band luminosity functions from z ~ 0-3, and slightly over- 
predict galaxies in the rest near-IR (K-band) at high redshift 
{z ~ 2-3). In all UV-NIR bands our models predict an excess 
of low-luminosity galax:ies , which confirms the excess found 
bv lFontanot et al.l l|2009al ) and others based on stellar mass 
function comparisons. 

It would have been unsurprising if our simple ap- 
proach for computing IR luminosities disagreed drasti- 
cally with more detailed radi ative transfer c a lculatio ns or 
with observations. However, iFontanot et al.l (|2009bl ) and 
iFontanot fc Somerville! (|2010h showed that implementing a 
recipe similar to the one we adopt here within the MOR- 
GANA SAM produced very similar results to the full ra- 
diative t ransfer calculations using the GRASIL code of 
ISilva et"a l. (1998), at least for global quantities such as lu- 
minosity functions and counts. And we have shown that 
our models reproduce galaxy counts from the UV to the 
mid-IR (-^ 70 fxm) remarkably well. A growing discrepancy 
starts to emerge at longer wavelengths: our models under- 
produce intermediate luminosity (~ 30-100 mjy) sources in 
the SPIRE 250 /im band by a factor of 2-5, and S ^ 5 mJy 
sources at 850 fim by an order of magnitude or more. This 
problem is far fr om being unique to our models, as we have 
discussed above. IClements et al.l (2010) show that none of 
the models that they compare with their SPIRE 250, 350 
and 500 /im count data agree with their results very well. 
Most of these models are empirical "backwards evolution" 
m odels, but th e y also compare with the semi-analytic model 
of iLacev et al.l ( 20081 ) . which predicts a significant excess of 
luminous galaxies in all three SPIRE bands. 

When comparing with estimates of luminosity functions 
at 2 ~ 0.5-2, we find deficits of luminous galaxies at high 
redshifts in the rest 8 /^m and 24 /im bands. In the mid-IR 
(8 and 24 /im), it is possible that there could be significant 
contamination from obscured AGN, particularly at z ~ 2 



(e.g. iDaddi et al.ll2007l ). The agreement with the estimated 
total IR LF to z ~ 2 is not terrible (within the observational 
errors). It is important to remember that these observational 
estimates rely on k-correcting from an observed wavelength 
to the rest-frame in a messy part of the SED (particularly 
in the case of 8 and 24 /im), or on estimating a total IR 
luminosity from a single or a limited number of observed 
wavelengths. These conversions are themselves highly un- 
certain and in general rely on SED templates. Therefore, it 
is encouraging that the agreement between our models and 
the more directly observed quantity, the galaxy counts, is in 
general superior to the agreement with the derived quanti- 
ties (rest-frame or total luminosity functions). In our com- 
panion paper (GSPD), we also show the redshift dependence 
of the build-up of the EBL at 24 /im, 70 /im, and 160 /im 
compared with available observations and find fairly good 
agreement. 

By integrating over all galaxies, accounting for the red- 
shifting and dilution of light, we estimate the integrated Ex- 
tragalactic Background Light predicted by our models. As 
expected, our EBL predictions lie close to the lower lim- 
its from integrated galaxy counts. The largest uncertainties 
(model-to-model differences) in our predictions are in the 
mid-IR, due to the limitations of available IR templates. 
These will improve as more data from multi-wavelength 
observations are synthesized. In particular, important con- 
straints on this part of the EBL, and correspondingly on 
the star formation history and dust SEDs of galaxies, may 
be obtained from observations of GeV and TeV gamma 
rays. High energy gamma rays are attenuated via electron- 
positron pair production against the EBL. In principle, the 
cosmological history of the EBL could be reconstructed by 
comparing observations of high-energy sources at different 
redshifts to their known intrinsic spectra. In a companion 
paper (GSPD), we provide a detailed analysis of the impli- 
cations of our predictions for current and future gamma ray 
observations. 
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